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FAST: A Fully Asynchronous Split Time-Integrator for Self- Gravitating 
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Abstract 

We describe a new algorithm for the integration of self-gravitating fluid systems using SPH method. 
We split the Hamiltonian of a self-gravitating fluid system to the gravitational potential and others (kinetic 
and internal energies) and use different time-steps for their integrations. The time integration is done in 
the way similar to that used in the mixed variable or multiple stepsize symplectic schemes. We performed 
three test calculations. One was the spherical collapse and the other was an explosion. We also performed 
a realistic test, in which the initial model was taken from a simulation of merging galaxies. In all test 
calculations, we found that the number of time-steps for gravitational interaction were reduced by nearly 
an order of magnitude when we adopted our integration method. In the case of the realistic test, in 
which the dark matter potential dominates the total system, the total calculation time was significantly 
reduced. Simulation results were almost the same with those of simulations with the ordinary individual 
time-step method. Our new method achieves good performance without sacrificing the accuracy of the 
time integration. 
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1. Introduction 

The number of particles used in simulations of 
galaxy formation with A^-body /Smoothed Particle 
Hydrodynamics (SPH) method has not increased much 
since the early days of Katz & Gunn (1991) and Navarro 
& Benz (1991), though the number of particles used 
in pure A^-body cosmological simulations has increased 
drastically. For A'^-body simulations, the largest run in 
1991 used - 2 X 10^ particles (Suto & Suginohara 1991) 
and the largest run recently performed used ^ 7 x 10^" 
particles (Kim et al. 2009). The number of particles 
has grown by nearly four orders of magnitudes in two 
decades. On the other hand, for A^-body/SPH simu- 
lations of galaxy formation, the first simulations used 
■~ 4000 SPH particles for a single halo (Katz & Gunn 
1991) and the largest simulation which is performed 
recently used ^ 3.2 x 10^ SPH particles for a single halo 
(Governato et al. 2009). ^ The scale up factor is only 80 
in two decades. This is because time-steps become quite 
short in dense and compact self-gravitating gas clouds of 
star-forming regions. 

This problem is severer in simulations with higher reso- 
lution, since these simulations resolve denser gas. In gen- 



eral, supernova (SN) explosion in dense regions leads the 
shortest time-step. Here we roughly estimate the decrease 
of the time-steps due to SNe. We consider a compact re- 
gion of the interstellar medium (ISM) with the tempera- 
ture of TisM, where the sound speed is cism, as a potential 
site of the star formation and that the region is rapidly 
heated to TgN, where the sound speed is Csnj by SN with 
the energy of i^sN- The contraction factor between the 
time-step of the ISM after the SN, dtsN, and before the 
SN, dtisM, in the compact region is 



dt^^/dti^yi — cism/csn, 

oc(Tism/Tsn)'/', 



(1) 



Note that a part of SPH particles were converted into star par- 
ticles, thus the number of SPH particles was reduced during the 
galaxy evolution. 



where i^sN is the energy of the single SN and m is 
the mass of the heated region or the mass resolution in 
Lagrange schemes such as SPH, respectively, and we use 
TsN Esi^/m. From this equation, we can easily find that 
the contraction factor becomes smaller when (i) mass reso- 
lution becomes higher and (ii) the temperature of the ISM 
becomes lower (see also section 2 for more detailed discus- 
sion). Thus high-resolution simulations which model the 
ISM with low temperature (< 10^ K) require much shorter 
time-steps than conventional simulations of galaxy forma- 
tion with a cooling cut off at ~ 10^ K. 

The individual time-step method (Aarscth 1963; 
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McMillan 1986; Makino 1991a) reduces the total calcu- 
lation cost significantly in simulations which cover a wide 
range of timcscalcs, by assigning different time-steps to 
different particles and integrating only a small fraction of 
particles with small time-steps. Here, we extend this idea 
for the time integration of self-gravitating fluid particles in 
order to achieve a further reduction of the total calculation 
cost. Our new method allows an individual fluid particle 
to have different time-steps for gravitational and hydro- 
dynamical interactions and integrates these interactions 
asynchronously. As stated earlier, the smallest time-steps 
are associated with particles heated by SNe feedback. 
These particles have the thermal and kinetic energy many 
orders of magnitudes larger than the gravitational poten- 
tial energy. Therefore, if we assign different time-steps 
to gravitational and hydrodynamical forces, we should be 
able to use much longer time-step for gravity, thereby ac- 
celerating simulations by a large factor. We named this 
time-integration scheme for self-gravitating fluid as FAST 
(Fully Asynchronous Split Time- integrator). 

There are two main advantages of the FAST method 
over the traditional individual time-step method for self- 
gravitating fluid simulations. First, FAST reduces unnec- 
essary gravitational force evaluations in small time-steps 
induced by SNe. Since the number of dark matter and 
stellar particles is usually larger than that of SPH particles 
in typical simulations of galaxy formation, the calculation 
cost of gravity is larger than that of hydrodynamics. This 
reduction of unnecessary evaluation of gravity is quite 
efficient for the acceleration of simulations. Simulations 
with hardware accelerators, such as GRAPEs (Sugimoto 
et al. 1990; Ito et al. 1991; Okumura et al. 1993; Makino 
ct al. 1997; Kawai et al. 2000; Makino et al. 2003), receive 
further benefit from FAST, since hardware accelerators 
are inefiicient in calculations with small number of par- 
ticles, because of small bandwidth and large latency of 
the bus between the host computer and the accelerator. 
The second advantage appears when we combine the in- 
dividual time-steps and the tree method (Barnes & Hut 
1986). Since the cost of the tree construction is indepen- 
dent of the number of active particles, it dominates the 
total calculation cost when the number of particles with 
small time-steps is small. Consequently, the total perfor- 
mance of simulation is not much improved by the use of 
individual time-steps. Table 1 of Wadsley et al. (2004) 
showed such a bad case. We can see that a half of the 
cost of smallest steps is that of the "Tree building" part. 
By adopting FAST, the number of tree construction is re- 
duced and hence simulations with the individual time-step 
and tree methods are accelerated significantly. 

There are several other ways to reduce the cost of sim- 
ulations with the individual time-steps and tree methods. 
McMillan & Aarseth (1993) applied the local update to 
the tree structure around particles which were updated, 
instead of reconstructing the whole tree structure at each 
step. This technique was used in GADGET-1/2 (Springel 
et al. 2001; Springel 2005). In VINE (Wetzstein et al. 
2009; Nelson et al. 2009), the construction frequency of 
tree structure was reduced by skipping several continuous 



time-steps and reusing old tree structure for force calcu- 
lation. They updated the tree structure at every ~ 10 
steps for the problem they showed in their paper. FAST 
method can be combined with these schemes to further 
reduce the cost of tree construction, if necessary. 

Our approach is similar to the multiple time-step 
method used in molecular dynamics, in which the long- 
range Coulomb force is updated less frequently than short- 
range van del Waals force (Streett et al. 1978). The main 
difference is that we combine the force splitting with in- 
dividual time-steps. 

The structure of this paper is as follows. In section 2, 
we estimate and compare the time-steps of particles in 
the hot region of star-forming galaxies. In section 3, we 
describe our new integration method for self-gravitating 
fluid, FAST. We briefly explain its implementation in §4. 
We present the results of test calculations and timing re- 
sults in section 5. A discussion on the maximum acceler- 
ation factor by FAST appears in section 6. In section 7, 
we provide summary. 

2. Estimate of Time-steps in Heated Regions of 
Star Forming Galaxies 

In this section, we estimate typical time-steps of an SPH 
particle in star forming regions of actively star forming 
galaxies in A^-body/SPH simulations of galaxy formation. 
This estimation allows us to estimate the maximum gain 
in the performance due to the use of the FAST scheme. We 
compare the typical Courant time-step of an SPH particle 
heated by SNe with the typical gravitational time-step of 
the particle. 

Here we estimate the typical Courant time-step of an 
SPH particle heated by SNe. For simplicity, we adopt fol- 
lowing four assumptions. First, we adopt a single stellar 
population (SSP) approximation for a star particle with 
Salpeter initial mass function (IMF) (Salpeter 1955) and 
the range of this IMF is set to be 0.1 Mq to 100 Af©. For 
this IMF, the speciflc SN rate is esN ^ 0.0074 SN/M©, 
where we assume 8 Mq or heavier stars become SNe at 
the final phase of their evolutions. Second, we assume that 
each SN injects the thermal energy of i?sN = 10^^ ergs to 
the surrounding ISM (the nearest iVNB particles). Third, 
we assume that the whole energy from SNe in a star par- 
ticle discharges in a single event (this is one of SN feed- 
back implementations proposed by Okamoto et al. 2008). 
Finally, we assume that the masses of the stellar and gas 
particles are the same. 

The mean additional internal energy for A'nb SPH par- 
ticles due to SNe of a single stellar particle is given by 

A'nb"t-sph 

~ 0.0074 X 10^1 X [ergs Mq-^, 

A'nb'tisph 

3.7 xlO^^ 1, 

[ergsg-i], (2) 

where to* and tosph are the masses of stellar and gas 
particles, respectively, and we use the relation m* =tosph- 



FAST 



3 



The sound speed, csn, of the heated gas region is 



CSN ^ 1)C^SN, 

6.4 X 10^ 



[km 



(3) 



where we assume an ideal gas with the adiabatic index of 
7 = 5/3. The original internal energy of the ISM before 
SNe, C/isM, is quite small, therefore we neglected ?7ism in 
the estimation of csn- The corresponding temperature of 
the heated region is TgN ~ 3.2 x 10^ [K]. 
Note that this temperature implies very short cooling 
timescale of ^ IQ^ yr. In real star-forming region, initially 
the SN ejecta have much high temperature, and the cool- 
ing time is much longer. In order to model SN feedback in 
a physically correct way, therefore, some tricks which pre- 
vent the quick radiative cooling (Gerritsen 1997; Thacker 
& Couchman 2000; Stinson et al. 2006) is necessary. We 
here assume some of these tricks are used. The size of an 
SPH particle, A, is 

3 mspH\^/^ 



A: 



47r 



(4) 



where p is the density of the SPH particle. Combining 
equations (3) and (4), we obtain the sound crossing time, 
tgN = A/cgN, in the heated region as follows: 



tsN ^ 4 X 10^ 



"T-SPH 



1/3/100 



1/3 



M, (5) 



. 1000 Mq . 

where we adopt Anb = 32 and Ah is the hydrogen number 
density of the heated region. The typical Courant time- 
step in the region is dfgN 0.1 x tgN- This equation tells 
us that the smallest time-step in simulations involving the 
low temperature ISM and SNe becomes shorter when mass 
resolution becomes higher and injected region becomes 
denser. 

By comparing equation (3) with the typical velocity of 
the ambient ISM, we can obtain the contraction factor of 
the time-steps caused by a SN explosion. Although the 
typical temperature of giant molecular clouds (GMCs) is 
low 10 K) and the corresponding sound speed in GMCs 
is small 0.2 km s^^), the linewidth of GMCs is higher 
than that expected by the sound speed of the ISM and 
predicts that GMCs are supported by supersonic turbu- 
lence. Thus we use empirical relations for the estimate 
of the timescale, instead of the local sound speed. The 
linewidth-size relation, which is often referred as Larson's 
law (Larson 1981; Solomon et al. 1987; Heyer & Brunt 
2004), gives us the typical velocity at the size of cloud. 
Larson's law is as follows: 

^cc.(^)^^^kms-], (6) 

where CTc and Lc are the linewidth and size of a cloud, 
respectively, and the applicable range of this relation is 
0.1 pc < Lc < 100 pc. Combining the virial theorem 
and this relation, we obtain cloud mass-linewidth relation 
(Solomon et al. 1987) that 



Afc = 2000 



V 1 km s" 



M, 



01 



(7) 



where Mc is a cloud mass. When we substitute A^NBTisPH 
into Mr, we obtain the mass resolution-linewidth relation: 



cTr = I - , 1 [km s J 



V 2000 Mq I ' 

This equation leads the velocity of the smallest cloud 
which can be expressed with the resolution of the sim- 
ulation. The contraction factor of the time-step in the 
ISM by a SN explosion, /cont, is 



CSN 



/ArNBTOSPH\l/'*/6.4x 102\-i 



V 2000 Mr. 



A^NB 



1/2 



l.SxlO^^f "^^P" 
VlOOO Mq 



1/4 



(9) 



where we adopted Anb = 32. This equation clearly shows 
that the Courant condition becomes quite tight in ISM 
heated by a SN explosion. When we use 6 km s~^, which 
is the sound speed of the ISM at lO'' K, as the typical 
velocity of the ISM, the contraction factor is /cont ~ 5.3 x 
10~^. We again adopted A^NB = 32. The contraction factor 
for simulations with the multiphase ISM and turbulence 
motions is smaller than that in conventional simulations 
of galaxy formation with a cooling cut off at 10^ K. 

In conventional simulations of galaxy formation, where 
the typical mass resolution is 10® Mq and the highest 
density of the ISM is 0.1 cm""^, the typical time-step for 
the heated region is dtsN 4 x 10^ yr. The typical grav- 
itational time-step, one-tenth of the local free-fall time 
at 0.1 cm~'^, is 5 X 10® yr. The difference between the 
Courant and gravity time-steps is ^ 10. On the other 
hand, in state-of-the-art simulations involving the multi- 
phase ISM, where tosph = 10^ Mq and Ar = 100 cm'^, 
the typical time-step is dtsN ~ 4 x 10'^ yr, whereas the the 
typical gravitational time-step at 100 cm"-^ is 1.6 x 10^ yr. 
The difference between two time-steps is ~ 40 and this 
difference is larger than that in conventional simulations. 
These simple estimates tell us that FAST can reduce grav- 
ity steps by a factor of 10 — 40. FAST is more efficient in 
simulations with high resolution. 

Thanks to the rapid increase of the computational 
power and the advance of numerical techniques, the mass 
resolution in current high resolution simulations of the 
galactic scale ISM has been quite high ('^ 1000 Mq). It 
will be soon reach the point than the mass resolution at 
where the number of SN events in an SSP particle is less 
than unity. ^ In such simulations, SN events are neces- 
sary to be treated as not an association of SNe in every 
SSP particle but discrete events in a fraction of SSP parti- 
cles so that the global SN event rate is consistent with the 
adopted IMF. Otherwise, we would introduce "fractional" 
SNe, which clearly would give wrong results for SNe feed- 
back. By this modification in the treatment of SNe, the 



When an SSP particle mass is lower than a critical mass mc, 
which is obtained by esN^c = 1, the number of SNe events in 
an SSP particle is lower than unity. If we use esN = 0.0074, 
rric ~ 135 Mq . 
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sound speed of the ISM in heated regions becomes much 
higher than that in equation (3) and the crossing time in 
these regions becomes much shorter than that in equation 
(5) . We show here a simple estimate of time-steps in the 
case where SN explosions are discrete events in SSP par- 
ticles. When we consider a SN explosion takes place in 
discrete manner, the received energy of the surrounding 
ISM of a compact region is modified as follows: 



10^1 

A^NB"lSPH 

5.0 X IQi^^l Mq 



ergs M0-1], 



, , [ergs g-^], (10) 

where we again neglected the original internal energy be- 
cause the value is sufficiently low compared with this 
value. The sound speed of the hot region is 



CSN,d = \Jl{l - l)C/sN,d, 

_ 7.5 X 10^/1 M0\i/2 
The contraction factor is 



(11) 



' A^NBmsPH \ r 7.5 X 103 / 1 Mq \ V 



V 2000 Mq ) L ATj^B^ 

: 2.7x10""^' > 



1/2 V77ISPH/ 



VIM,; ■ <>^) 

Note that the mass dependency in this equation is much 
stronger than that in equation (9). Combining equations 
(11) and (4), the sound crossing time in the heated region, 
tsN.d, is 



tsN,d^3.3x 10 



2/"^sph\^/^/100 cm 3n^i/3 



M, (13) 



where we adopt A^nb — 32. Wc find that the mass res- 
olution dependence in the equation (13) is stronger than 
that in the equation (5). Thus high resolution simulations 
of near future will be much harder than those of present. 
For efficient simulations, we have to introduce efficient nu- 
merical techniques which can handle a very wide range of 
time-steps. We believe that our new scheme will play an 
important role not only in current simulations but also in 
new simulations of galaxy formation in the near future. 

3. Basic Idea 

The basic idea of our new scheme is as follows. We 
allow gas (SPH) particles to have different time-steps for 
gravitational and hydrodynamical integrations. Thus, we 
extend the idea of individual time-steps, which allows dif- 
ferent particles to have different time-steps, to allow sin- 
gle particle to have different time-steps for different in- 
teractions. We then asynchronously integrate gravity and 
hydrodynamics with these different time-steps. This is 
the essence of our FAST method. Since the problem we 



have to solve is the time-integration of very hot gas par- 
ticles formed by SNe, we allow time-steps for gravity to 
be longer than those for hydrodynamics. To allow differ- 
ent time-steps for gravity and hydrodynamics, we use the 
technique of constructing multi-timestep symplectic inte- 
grator. We divide the Hamiltonian of a self-gravitating 
fluid system into a gravitational potential term and oth- 
ers, and integrate each part with its own time-step. 

The Hamiltonian of a self-gravitating fluid system of N 
gas particles is expressed as 



N 

E 



Vi 
2mj 



s) - 



N „ 



(14) 



where pi and qi are conjugate variables of the canoni- 
cal equation for particle i, is the mass of particle i, 
U is the internal energy of fluid, which is a function of 
<7, density, p, and entropy, s. Here, q, p, and s denote 

(91,<?2,93,---,9JV), {pi,P2,P3,---,Pn), and (si,S2,S3,---,SAr), 

respectively. Since we take into account arbitrary forms 
of hydrodynamical interactions, we express the internal 
energy for fluid as U{q,p,s). The first, second, and third 
terms in the right hand side of equation (14) are the ki- 
netic, internal, and gravitational potential energy of the 
system, respectively. The actual equations for p and s 
contain the contributions of non-conservative terms like 
artificial viscosity and radiative cooling/heating. For sim- 
plicity, we here regard the system as adiabatic (i.e., Si 
are treated as constants). Hence the internal energy term 
becomes the function of (q,p) and can be regarded as a 
potential term in the Hamiltonian. 

We split the Hamiltonian into the gravitational poten- 
tial term and others (see appendix 1): 



N 



H, 



grav 



N 

■V ^ Grriimj 



(15) 



(16) 



Wc then obtain the following expression of a symplectic 
integrator with the second-order accuracy, 

/(t + Af) «e^i'«-->e^*{-«--°>e^i'«-->/(t), (17) 

where "{ , }" is a Poisson bracket and At is a time-step. 
The equation (17) can schematically rewrite as follows: 



Vq = Vo + -At Ograv, 

x'o — > (Hydro update) 
v'q — > (Hydro update) 



vi = v[ + - At 



''grav 1 



(18) 

(19) 
(20) 

(21) 



where x, v, v' , and Cgrav indicate the position, the veloc- 
ity, the half-step advanced velocity, and the acceleration 
of gravitational force, respectively. Subscripts and 1 
indicate epochs of time-integration at t and t + Ai, re- 
spectively. 
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There are many ways to integrate the hydrodynamical 
part of equation (17), we here choose the second-order 
symplcctic method (e.g., Hernquist & Katz 1989). We 
divide equation (15) again into the foUowing two parts: 

N 2 
-ffhydro^T 2^' 

Hhydro,u ^U{q,p). (23) 

Consider the case that Atg = lAth, where I is a natural 
number. We obtain a new expression of equation (17) as 

/(t + At)«e^{'^'=-> 

Jg — 2^ { I -Whydro , U } g A* (i { , hydro , T } g — 2^ { ^ -f^hydro , U } J ' 

e^^'^'=->/(t). (24) 

This equation tehs us that we can reduce the computa- 
tional cost of gravity if Atg > Ath {I > 1). If we adopt 
1 = 1, the integrator is the same as the standard "leap- 
frog" method for self-gravitating fluid. 

In figure 1, we show schematic pictures of the usual 
leap-frog and FAST methods. For FAST, we consider the 
case that digrav = 2 dihydro- The computational cost of 
gravitational force in FAST becomes half of the leap-frog 
method in this case. In practice, the time-step ratio, I, 
adaptively changes. 

It should be noted that, even though wc borrowed the 
formalism of symplcctic integrators to describe our FAST 
method, the FAST method itself is not symplcctic. This is 
because we change the time-steps for gravitational and hy- 
drodynamical interactions, after we split the Hamiltonian. 
In addition, we use different time-steps for different par- 
ticles. However, this issue is not as crucial as that is for 
pure A^-body simulations, since the usual hydrodynami- 
cal simulations introduce a dissipation term. The time- 
integration of hydrodynamic simulations is usually time 
irreversible and is breaking the symplcctic nature inher- 
ently. 

There have been several proposed methods which can 
retain either symplecticness (Farr & Bertschinger 2007) or 
time symmetry (Makino et al. 2006) when used with the 
individual time-step method. However, these schemes are 
computationally expensive and it is not clear if the use 
of these schemes is worthwhile or not. In this paper, we 
concentrate on the traditional, non-symplectic implemen- 
tation of individual time-step algorithm and its extension. 

4. Implementation 

4-.1. The Code 

The code used in this paper is a parallel tree SPH 
code, ASURA, which utilizes the special purpose hardware 
GRAPE (Saitoh in prep.). Gravitational force was solved 
by Tree with GRAPE (Makino 1991b). In this paper, 
we used the Phantom GRAPE library for calculations of 
gravity, which is a software emulator of GRAPE pipelines 
(kindly provided by Kohji Yoshikawa). We used an open- 
ing angle of 0.5 and only monopole moments for force 



calculations. Hydrodynamics was followed by the stan- 
dard SPH method (e.g., Lucy 1977; Gingold & Monaghan 
1977; Monaghan 1992). Wc used the "gather" formulation 
of SPH for the density estimation, whereas the "gather 
and scatter" formulation of SPH for the pressure gradi- 
ent and the time derivation of internal energy (Monaghan 
1992). We adopted the asymmetric form energy equation 
(e.g., Steinmetz & Mueller 1993). We iteratively deter- 
mined the kernel radius of each SPH particle in every step 
in order to keep the number of neighbor particles, 32 ± 2. 
We used an artificial viscosity term, of which form is the 
same as that proposed by Monaghan (1997), in order to 
handle shocks. The value of the viscosity parameter, a, 
was set to be unity. ASURA adopts the variable and in- 
dividual time-step method (McMillan 1986; Hernquist & 
Katz 1989). Following Makino (1991a), ASURA adopts an 
extended version of the individual time-step method, i.e., 
the "hierarchical" time-step method where time-steps are 
quantized by the power of two of a baseline time-step in 
order to improve the simulation performance with individ- 
ual time-steps. We also implemented the time-step limiter 
for hydrodynamical interactions in order to maintain the 
difference of time-steps among neighbor particles small 
enough (Saitoh & Makino 2009). Here we adopted the 
factor of the time-step difference in neighbors, / = 4. The 
current version of ASURA implements two time integrators, 
namely the ordinary leap-frog and FAST methods. 

4.2. Time-steps 

The time-steps were determined as follows. The gravity 
time-step of an i-th particle was estimated by 

rftgrav.i = C'grav minf , /^j^^^ -, ] "'^'''"'^ ' ' , (25) 

where e is a gravitational softening length, Cgrav is a pa- 
rameter which controls the accuracy (we here adopt 0.1), 
and agrav,j is the time derivation of the acceleration, re- 
spectively. 

Following Monaghan (1997), the hydrodynamical time- 
step of the i-th SPH particle was determined by 

rf^hydro,j = C'hydro — , (26) 

^sig.i 

where is the kernel size of the SPH particle (the inter- 
action scale is 2hi), Chydro = 0.25, and Wsig,i is the local 
maximum signal- velocity of i-th particle defined by 

Vsig,i = max(ci + Cj - 3w,j ) , (27) 

3 

where j indicates the indices of neighbor particles, Ci and 
Cj is the sound speed of i-th and j-th SPH particles and 
Wij = Vij ■ Xij/\xij I is a projected relative velocity between 
the SPH particles. We set Wij = if Wij > 0. 

When we use the FAST method, we asynchronously 
integrate gravity and hydrodynamics by the leap-frog 
method with different time-steps for gravity (Eq. 25) 
and hydrodynamics (Eq. 26). If digiav 7^ 2"o?thydro where 
n is an integer number of > 0, we change the gravita- 
tional time-step so that it satisfies the above criterion 
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Leap-frog 



Drift (free motion) , Drift (free motion) , Drift (free motion) , Drift (free motion) 



dt 




dt 




dt 




dt 



Kick Kick Kick Kick Kick 

(gravity+hydro) (gravity+hydro) (gravity+hydro) (gravity+hydro) (gravity+hydro) 



FAST (dtgrav = 2 dthydro) 



Drift (free motion) Drift (free motion) , Drift (free motion) Drift (free motion) 




dthydro .^^^ dthydro 

Kick(hydro) 



(gravity+hydro) 



dt, 



grav 




dthydro -^^^ dthydro 

Kick(hydro) 



(gravity+hydro) 



dt, 



grav 




(gravity+hydro) 



Fig. 1. The schematic picture of the leap-frog and FAST methods for the integration of a self-gravitating 
fluid. See also figure 1 in Fujii et al. (2007) for MVS and BRIDGE methods. "Kick" means the mo- 
mentum exchanges between particles, while "Drift" denotes the free (inertial) motions under given velocity vectors. 



in the following way: we reduce the time-step of grav- 
ity to digj-av, where dtgrav > f^^grav = 2"c?thydro, and n is 
the maximum integer number that satisfies this relation. 
If dtgrav < dthydro- wc rcducc dthydro to the same value as 

dtgrav ■ 

When we used the ordinary leap-frog method, we picked 
up the smaller one of the two time-steps as the time-step 
of an SPH particle. 



dt = min(dtgrav, dthydro), 



(28) 



and synchronously integrated both gravity and hydrody- 
namics. In general, the acceleration and its differential 
terms in equation (25) should be measured relative to the 
total acceleration {i.e., the sum of the gravitational and 
hydrodynamical accelerations). However, in the hydro- 
dynamical simulations, the Courant condition leads the 
smaller time-step compared with the time-step obtained 
by the total acceleration. Therefore the simple determi- 
nation by equation (28) worked sufficiently. 

5. Numerical tests 

We performed three tests. The first test was the col- 
lapse of a gas cloud and the second test was the point-like 
explosion of a self-gravitating gas cloud. The third test 
was a more realistic simulation. We performed simulations 
of galaxy-galaxy collisions, where galaxies consist of dark 
matter, star and gas particles. These tests incorporated 
both gravity and hydrodynamics and were representative 
of the evolution of self-gravitating fluid in galaxy forma- 
tion simulation or other astrophysical simulations. We, 
hereafter, denote the results with the ordinary individual 
time-step method as "Ind" and the results with the indi- 
vidual time-step with the FAST method as "FAST", in 
tables and figures. The first two tests were designed as 
simple tests for the validity of the FAST method, while 



in the third test we investigated the actual gain in the 
calculation speed as well as the accuracy of the result. 

The first and second tests were done on a system with 
a 2.4 GHz Opteron 280 processor (Italy core), while the 
third test was done on 2.2 GHz quad-core Opteron pro- 
cessors (Barcelona core) of Cray XT4 system at Center 
for Computational Astrophysics of National Astronomical 
Observatory of Japan. We used one CPU core for the first 
and second tests, whereas we used 128 CPU cores for the 
third test. 

5.1. Test I: Three dimensional self- gravitational collapse 
tests 

We performed the integration of three-dimensional 
spherical coUapse of adiabatic gas (e.g., Evrard 1988; 
Hernquist & Katz 1989). This test is one of standard 
tests for SPH method which involves self-gravity. 

We prepared a gas sphere with the total mass and the 
radius both unity. The gravitational constant was also set 
to be unity. The initial profile of the gas sphere was p(r) oc 
1/r, where r is the distance from the center of coordinates. 
The adiabatic index and the specific internal energy of the 
gas were set to be 7 = 5/3 and 0.05, respectively. The 
gas sphere had a negative value of the total energy, E ~ 
—0.6. When the evolution starts, the sphere begins to 
collapse. The shock takes place in the central region and it 
propagates outward. Finally, the system reaches the state 
of virial equilibrium. In this test, we used 30976 particles 
for the sphere and we followed the evolution of the gas 
sphere to T = 3. We adopted 0.038 for the gravitational 
softening length. 

Figure 2 shows radial profiles of density, pressure, and 
radial velocity in three different epochs, T = 0.9, 1.2 and 
2.4, for both the ordinary individual time-step and FAST 
methods. In this figure, we plotted mean physical quanti- 
ties of every 300 particles. It is obvious that the results of 
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two methods are identical. We also confirmed that these 
results agree well with the result obtained using global 
time-steps. Therefore we can conclude our new method is 
accurate enough. 

Figure 3 shows the values of time-steps for gravity and 
hydrodynamics for the run with the FAST method as 
a function of the distance from the center at T = 0.9. 
Time-steps for gravity and hydrodynamics are different 
in the post-shock region and the same in the ambient, 
pre-shocked region. The transient region clearly matches 
with the shock front at the radius of ^ 0.2 (see the top- 
left panel of figure 2). The values of hydro and dtg^-^^ in 
the post-shock region differ by a factor up to four. Figure 
4 shows cumulative fractions of dthydro and dtgrav for the 
simulation with the individual time-step with FAST. We 
can see that almost half of the particles have hydro time- 
steps smaller than the minimum time-step for the gravity. 
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Fig. 3. Radial profiles of dthydro and ^igrav for the 
spherical collapse test with the individual time-step with 
FAST. The epoch is T = 0.9. Sohd and dotted 
curves indicate dfgrav and dthydroi respectively. The 
hatched region corresponds with the soften region for 
the particle at the center by the gravitational softening. 



The errors of the total energy, the difference between 
the values of the total energy at the initial (T = 0) and 
final (T ~ 3) states, for the Ind and FAST methods are 
shown in table 1. Eg and Ef represent the initial (T = 0) 
total energy and the final (T = 3) total energy, respec- 
tively. These values are acceptably small. In our test 
runs, the absolute value of the energy error for the time 
integration with FAST is smaller than that for the time 
integration without FAST even though the time-step for 
gravity is larger. The errors caused by the gravitational 
and hydrodynamical integrations had the opposite signs 
and partially canceled each other. 

In table 2, we show timing results of the collapse test. 
Our new method is faster than the original method, but 
not by a large factor. The calculation time of gravity is 
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Fig. 4. Cumulative fractions of particles as a func- 
tion of rffhydro and dtgrav for the spherical collapse 
test with the individual time-step with FAST. The 
epoch is T = 0.9. Solid and dotted lines indicate 
cumulative fractions of dtgrnv and dthydroi respectively. 



Table 1. Total energy errors for the spherical collapse test. 



Method 


\iEs-Ef)/Es\ 


Ind 


3.0 X 10-^ 


FAST 


1.3 X 10-3 



reduced to two thirds, and that of tree construction is 
reduced to two fifths. However, in this test the hydrody- 
namics part dominates the total cost. 

Table 3 shows the number of steps and integrated par- 
ticles for the collapse test. The ordinary individual time- 
step method required 1979 steps for the simulation in our 
implementation. The FAST method required 2022 steps 
for the hydrodynamics part and 778 steps for the grav- 
ity part. The reduction of the calculation time for tree 
construction is directly proportional to the reduction of 
gravity steps. In our new method, the total number of 
integrated particles for the gravity part becomes almost 
two thirds of that for the hydrodynamics part in this test. 

5.2. Test II: Three-dimensional explosion tests 

Now, we discuss the result of three-dimensional explo- 
sion test. We designed this test to mimic an explosion of 

Table 2. Timing results for the spherical collapse test. 







Time [sec] 




Method 


Total 


Gravity" Hydro 


Others 


Ind 


1754 


682 (52) 960 


112 


FAST 


1523 


468 (20) 943 


112 



" The tree structure construction times in the gravity part 
are shown as parenthetical numbers. 
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Fig. 2. Radial profiles of density (left), pressure (mid), and radial velocity (right) at T = 0.9, 1.2, and 
2.4 (top to bottom). Horizontal axis is the distance from the origin of coordinates. Curves and cir- 

cles indicate the profiles obtained with the individual and synchronous time-steps for gravity and hydro- 
dynamics, "Ind", and the individual and asynchronous time-steps for gravity and hydrodynamics, "FAST" . 
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Table 3. Steps 
particles for 



and 
the 



number 
spherical 



of integrated 
collapse test. 





Gravity 


Hydro 


Method 


steps A^int.grav 


steps iVint, hydro 


Ind 


1979 2.4 X 10^ 


1979 2.4 X 10'^ 


FAST 


778 1.7 x 10^ 


2022 2.4 X 10^ 



a single SN in a self-gravitating gas cloud. When we take 
the mass of 10^ Mq and the radius of 100 pc as typical 
values of a giant molecular cloud (Dame et al. 1986), its 
potential energy is ~ 10^"'^ ergs (here we assume the cloud 
is in a virial equilibrium state). This value is comparable 
to the energy released by a single Type II SN. Therefore, 
in order to investigate the behavior of an exploding cloud 
induced by SN, we solved the evolution of a gas cloud with 
a positive total energy comparable to the absolute value 
of the original total energy. 

We used the particle distribution of the three- 
dimensional collapse test at T = 3 as the initial particle 
distribution of this explosion test. We added the thermal 
energy in the central 32 particles with SPH manner. Since 
the original total energy of the system is ^ — 0.6, the new 
total energy of the system was set to be £' = 0.5,1,2, and 
10. We set T to zero before the first step of the explosion 
calculation and follow the evolution to T = 5. In this test, 
we also plot the results of the global time-step case. 

Figure 5 shows snapshots of the expanding cloud, for 
the case oi E = 2, obtained with the FAST method at six 
different epochs (T = 0,1,2,3,4 and 5). The particles in 
the region \z\ < 0.1 are shown in this figure. The initially 
compact gas cloud expands driven by the high pressure 
gas added in the center of the cloud, and forms a spherical 
shell-like structure. The shell moves outward, and at the 
final phase (T = 5), the radius of the shell becomes ^ 5. 

Figure 6 shows time evolutions of density peaks for sim- 
ulations with several different values of the injection en- 
ergy. The positions of peaks are derived by averaging the 
positions of the 10 particles with highest local density. 
For the reference, in this figure, we plotted the results 
obtained with the global-step method. There are good 
agreements between individual time-steps with/without 
FAST runs and the global time-step runs. This is be- 
cause we adopted the time-step limiter for hydrodynamics 
(Saitoh & Makino 2009). Without this limiter, we would 
have failed to obtain agreements between different meth- 
ods. The difference of the the positions between individ- 
ual time-steps with/ without FAST runs are summarized 
in table 4. i?ind and i?FAST represent radii of shells at 
T = 5 for individual time-step without /with FAST and 
global time-step runs, respectively. In this table, we also 
show the difference between i?FAST and i?Giobai, which is 
the radius of the shell at T = 5 for the global time-step 
runs. The difference between FAST and Ind is compara- 
ble or smaller than the difference for the result of global 
time-step. 

Table 5 shows the timing results for the explosion test. 
We can see that the reduction in the cost of gravity calcu- 
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Fig. 6. Positions of density peaks as a function of time for 
various total energy cases. Solid, dashed, and dotted lines 
indicate evolutions of density peaks for cases of the global 
time-step (Global), ordinary individual time-steps (Ind), and 
individual time-steps with FAST (FAST). Numbers just 
above the lines indicate the values of the total energy. 



Table 5. Timing results for cloud explosion tests. 









Time 


[sec] 




Method 


E 


Total 


Gravity" 


Hydro 


Others 


Ind 


0.5 


1264 


498 (138) 


550 


216 


FAST 


0.5 


972 


175 (16) 


568 


229 


Ind 


1 


1034 


396 (100) 


458 


180 


FAST 


1 


799 


135 (15) 


482 


182 


Ind 


2 


941 


358 (92) 


438 


145 


FAST 


2 


719 


112 (14) 


451 


156 


Ind 


10 


1018 


383 (106) 


466 


169 


FAST 


10 


758 


89 (13) 


484 


185 



Items are the same as table 3. 

lation is much larger than that in the collapse test, and the 
reduction in the cost of tree construction is even larger. 
For E = 10, reduction in the tree construction cost is a 
factor of eight. 

The number of integrated particles and steps are sum- 
marized in table 6. By using the FAST method, we can 
greatly reduce steps for the gravity part. In these simula- 
tions, the use of the FAST method resulted in the reduc- 
tion of the number of gravity steps by a factor of 7 to 9. 
The numbers of integrated particles for the gravity part 
are reduced to only 0.5 — 0.7 times that of Ind runs. As 
is shown above, the speed up factors for the gravity part 
are 2.8 — 4.3. These results indicate that the decrease of 
the number of steps is quite efficient for integrations of 
self-gravitating fluid. There is almost no change in the 
hydrodynamics part. 
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Fig. 5. Snapshots of the expanding cloud at six different epochs (T = 0, 1, 2, 3, 4 and 5). Projected particle distribu- 
tions in a thin (1^1 < 0.1) region are shown. Dots indicate projected particle positions. This is the case that the ex- 
plosion simulation with E = 2. The time-integration was done by the individual time-steps with the FAST method. 



Table 4. Differences between the peak positions with individual time-step with FAST and others (individual time-step without FAST 
and global time-step). 







£; = 0.05 


E = 0.l 


E = 0.2 


£; = 10 


1 Rln 

1 -Rfast - 


-«:„.| 

--Rciobal 1 


1.5 % 
6.4 % 


1.6 % 
6.2 % 


1.3 % 
1.9 % 


0.4 % 
1.1 % 


1 Rn 





Table 6. Steps 


and 


number 




of inte- 


grated 


particles for 


cloud 


explosion tests. 






Gravity 




Hydro 


Method 


E 


steps 


A^int,grav 


steps 


^ int, hydro 


Ind 


0.5 


4675 


9.0 X 10^ 


4675 


9.0 X 10*' 


FAST 


0.5 


525 


5.9 X 10*^ 


4431 


8.8 X 10*^ 


Ind 


1 


3480 


7.8 X 10^ 


3480 


7.8 X 10" 


FAST 


1 


479 


4.8 X 10^ 


3535 


7.8 X 10" 


Ind 


2 


3515 


7.5 X 10« 


3515 


7.5 X 10" 


FAST 


2 


452 


4.1 X 10^ 


3397 


7.4 X 10" 


Ind 


10 


3104 


7.5 X 10" 


3104 


7.5 X 10" 


FAST 


10 


427 


3.4 X 10^ 


3412 


7.6 X 10" 



5.3. Test III: Merger simulations 

In this section, wc discuss the rcsuh of the apphca- 
tion of the FAST method to a reahstic problem, namely 
simulations of galaxy-galaxy collisions. Simulations we 
performed here were based on our recent galaxy-galaxy 



collision simulations of Saitoh et al. (2009), in which we 
followed the cooling of gas down to 10 K. We used the 
model of MIC. Gravity, hydrodynamics, radiative cooling, 
far-ultraviolet heating, star formation, and type-II SNe 
were taken into account. The condition for the star for- 
mation is that the gas is dense (tih > 100 cm^'^) and cold 
(T < 100 K) with converging flows. The regions which sat- 
isfy these conditions form stars following the Schmidt-law 
with the local star-formation efficiency of 0.033. Further 
details of the modeling of star formation were described in 
Saitoh ct al. (2008) and Saitoh ct al. (2009). Gravitational 
softening was set to be 20 pc for all particles. The initial 
numbers of dark matter, (old) star, and gas particles were 
6930000, 341896, and 148104, respectively. We used 128 
cores of Cray XT4 system at Center for Computational 
Astrophysics of National Astronomical Observatory of 
Japan. 

Figure 7 shows density and temperature maps for 
merger simulations at T = 420 Myr by individual time- 
steps without and with the FAST method. Wc can easily 



No. 



FAST 



11 



see that these two integration methods show quite sim- 
ilar resuhs in density and temperature structures. The 
positions of "Heat spots" due to SNc arc not perfectly 
identical because of run-to-run fluctuations. Other global 
properties of these galaxies are also identical for both runs. 
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Fig. 8. Cumulative fractions of particles as a function of 
'^thydro S'lid ^igrav for SPH particles and dt^^^dy for col- 
lisionless particles. The time-steps are sampled from the 
merger simulation by the individual time-step with FAST. 
The epoch is T = 0.43 Gyr. Solid, dotted, and dashed his- 
tograms indicate cumulative fractions of dtgrav and dthydro 
for SPH particles and dtnbody for coUisionlcss particles. 
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Fig. 7. Density and temperature maps for simulations by 
individual time-steps without and with the FAST method. 
Each panel shows 16 kpc X 16 kpc in the orbital plane. 
The left and right columns show the results by individ- 
ual time-steps without and with the FAST method, re- 
spectively. The epoch of these maps is T = 420 Myr. 

Figure 8 shows cumulative fractions of particles as a 
function of dihydro and rfigrav for SPH particles and dtnbody 
for coUisionlcss particles. The minimum time-step for the 
hydrodynamics part is shorter than that for the grav- 
ity part for SPH particles by a factor of eight. In ad- 
dition. coUisionlcss particles have longer time-steps than 
SPH particles. Therefore the gravity part is skipped in 
the lowest 3 levels and the calculation is accelerated sig- 
nificantly. 

Table 7 shows timing results of merging simulations. 
We sampled two typical epochs, i.e., 350 Myr < T < 
400 Myr and 400 Myr < T < 450 Myr. The former epoch 
is a quiescent star forming phase before the first encounter 
while the later epoch is a significantly enhanced star form- 
ing phase during the first encounter. The total integration 
time for the simulation with the FAST method decreases 
by almost a factor of two from that of the simulation with- 
out the FAST method. With the FAST method, the grav- 
ity part is ^ 7 times faster than that without the FAST 
method. The reduction of the total calculation time is 
similar for quiescent and starburst phases. 



In table 8, we show numbers of time-steps and inte- 
grated particles for gravity and hydrodynamics parts. The 
simulation with the FAST method required 7 times 
smaller number of gravity steps than that without the 
FAST method. Note that the number of integrated parti- 
cles for gravity is almost the same for the Ind and FAST 
method. The large reduction in the calculation time is due 
to both the reduction of the number of tree constructions 
and the removal of force calculations with small number 
of particles, where the calculation becomes inefficient in 
individual time-steps or parallel computers. 

6. Maximum Acceleration Factor by FAST 

In this section, we estimate the maximum acceleration 
factor due to the introduction of the FAST method using 
a simple calculation cost model. By comparing the cal- 
culation costs of the runs with/without FAST, we obtain 
the acceleration factor due to the FAST method. 

We here model the calculation cost of a simulation as 
follow: 

iall = ^s,g(<t,g + N^Jcs) + A^s,h(it,h + ^u,hto,h), (29) 
= iVs_g<gi.av + A^'s^h^hydro- (30) 

where Ns^g and Ng.h are the number of steps for grav- 
ity and hydrodynamics, tt,g and tt,h arc the calculation 
times of tree construction for gravity and hydrodynamics, 
A^u.g and Nu,h are the mean number of updated parti- 
cles for gravity and hydrodynamics in each step, <e,g and 
<c,h are the mean evaluation time of gravity and hydro- 
dynamical interactions for a single particle, respectively. 
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Table 7. Timing results for merger simulations. 









Time 


[sec] 




Method 


Epoch 


Total 


Gravity" 


Hydro 


Others" 


Ind 
FAST 


350 Myr 400 Myr 
350 Myr ^ 400 Myr 


14301 
7249 


5887 (2132) 
919 (264) 


3545 
2923 


4869 
3407 


Ind 
FAST 


400 Myr 450 Myr 
400 Myr -> 450 Myr 


16454 
9646 


6703 (2441) 
953 (279) 


4041 
4152 


5710 
4541 



" Items are the same as table 3. " In this test runs, "Others" includes the calculation times of the radiative cooling, 
star formation, and SNe routines. 



Table 8. Steps and number of integrated particles for merger simulations. 







Gravity 


Hydro 


Method 


Epoch 


steps 


-^int,grav 


steps 


^ int. hydro 


Ind 


350 Myr ^ 400 Myr 


18697 


1.6 X lO'J 


18697 


5.8 X 10'' 


FAST 


350 Myr 400 Myr 


2310 


1.6 X 10^ 


16626 


6.4 X 10^ 


Ind 


400 Myr 450 Myr 


21296 


1.8 X 10« 


21296 


7.3 X 10' 


FAST 


400 Myr ^ 450 Myr 


2425 


1.8 X 10^ 


22959 


8.0 X 10^ 



and tgrav and thydio arc the mean calculation times in a 
single step for gravity and hydrodynamics, respectively. 
In this model, we neglected the calculation cost of miscel- 
laneous operations, such as domain decomposition, time- 
integration, evaluations of time-steps. 

In traditional individual time-steps, the number of grav- 
ity steps is the same as that of hydrodynamical steps. 
Therefore the total calculation cost with traditional indi- 
vidual time-steps is 

tall, Ind = Nsh ,Ind (^grav.Ind ro.Iiid )• (31) 

On the other hand, in FAST, the number of gravity steps 
is different from that of hydrodynamical steps. According 
to the argument in section 2, the number of hydrody- 
namical steps is about ten times larger than the gravity 
steps. It means A's h^pAST = 10 x iVs g PAST in equation 
(30). Thus, the total calculation cost with FAST expresses 
as 

^all.FAST ~ -^s^h,FAST(Y^tgrav,FAST + ^hydro.FAST)- (32) 

The acceleration factor is defined as t=: tail. ind/^aii, fast- 
If we assume that the number of hydrodynamical steps in 
the calculation without FAST is the same as that with 
FAST and that calculation times between with/without 
FAST are the same, the acceleration factor is 

10(tgiav + ^hydro) 

'^ = 1 TTnT ■ (^^^ 

''grav ~r -L^f'hydro 

In this equation, we removed suffixes Ind and FAST of 
calculation costs. Here we consider three typical cases 

that tgrav ^hydro; ^grav — ^hydro; and tgrav ^hydro- The 

first case corresponds to usual 7V-body/SPH simulations. 
In this case, the acceleration factor is r = 10, and is quite 
large. The second case corresponds to the case that (a) the 
calculation time of the gravitational force reduces signif- 
icantly by adopting hardware/software accelerators, such 



as GRAPE or Phantom-GRAPE, and/or (b) the calcu- 
lation cost of hydrodynamics is rather expensive because 
of the treatment of complex baryon physics, for instance 
star formation, SNe, and chemical evolution. The accel- 
eration factor in this case is r 2. Our simulations are 
close to the second case. The final case is the ideal case 
that the calculation cost for gravity is negligible. In this 
case, the acceleration factor becomes unity and the gain 
due to FAST is zero. We do not think this hypothetical 
situation can occur. 

7. Summary 

In this paper, we describe a fast integrated method, 
"FAST" for self-gravitating fluid. The FAST method as- 
signs different time-steps for gravitational and hydrody- 
namical interactions and integrates them asynchronously. 
The formulation of the FAST method is similar to the 
multi time-step method (Streett et al. 1978) and also 
regarded as an extension of "multi-step" symplectic in- 
tegrators, such as mixed variable symplectic (Wisdom 
& Holman 1991), multiple stepsize (Skeel & Biesiadccki 
1994), and the BRIDGE (Fujh et al. 2007) methods. 

The approach of the FAST method is qualitatively dif- 
ferent from other reduction techniques of the tree con- 
struction in gravity part. Thus the FAST method elimi- 
nates unnecessary tree constructions and gravity calcula- 
tions by adopting longer time-steps for gravitational evo- 
lution than that for hydrodynamics. 

We found that the evolution of collapsing and exploding 
self-gravitating fluid by the FAST method are identical to 
these by the usual unsplit method which integrates gravity 
and hydrodynamics synchronously. 

As a realistic test, we applied the FAST method to 
merger simulations including self-gravity, hydrodynamics, 
radiative cooling, far-ultraviolet heating, and SN (Saitoh 
et al. 2009). In this test, we found that simulations with 
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and without the FAST method showed quite similar evo- 
lution. The calculation with FAST was nearly a factor 
of two faster. This large gain was due to the reduction 
in the gravity steps with small number of particles. The 
FAST method is very effective in accelerating simulations 
of self-gravitating fluid. 
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Appendix 1. 
its variants 



Symplectic integration method and 



In this appendix, we explain the symplectic integration 
method briefly. Then we explain sophisticated versions of 
symplectic integration methods, i.e., mixed variable sym- 
plectic method (Wisdom & Holman 1991), multiple step- 
size method (Skeel & Biesiadecki 1994), and the BRIDGE 
(Fujii et al. 2007) methods. 

Symplectic integration methods (e.g., Dragt & Finn 
1976; Forest & Ruth 1990; Yoshida 1990; Yoshida 1993) 
are now widely used in the simulations of gravitating A'^- 
body systems. These methods preserve the symplectic 
form of the canonical equation of motions when calculat- 
ing the time variation of the system. This character leads 
to the very good conservation of system's total energy on 
the course of numerical calculation. 

When we express H as the Hamiltonian of the system, 
and p and q as six-dimension coordinates, canonical equa- 
tions are 



dq 
dt 
dp 
'di 



dH 
dp ' 
dH 
dq 



We can summarize above two equations as 



dl_ 
dt 



(Al) 
(A2) 

(A3) 



where / is p or 5, respectively, and {,} is a Poisson bracket. 
We define an operator that 



lH}f. 



(A5) 



We can write a generalized canonical equation (A3) as 
dt 

When we integrate equation (A5) from t to t + At, the 
formal solution of the equation (A5) is written as 

/(t + At) = e^'{'^>/(t). (A6) 

Here, we consider the Hamiltonian of a self-gravitating 
system with particles. In this case, the Hamiltonian is 
written as 



H = Ha + Hb, 
where 

N n 



N 



Grriim j 



qij 



The formal solution is written as 
/(t + Ai) = e^*«'«*}+{^^-»/(i). 



(A7) 

(A8) 
(A9) 

(AlO) 



Applying Barker-Champbell-Hausdorff formula 
(Varadarajan 1984) to equation (AlO), we obtain a 
first order integrator 



/(t + Ai)«e^'{'^*>e^*^'^«>/(0, 
and a second order integrator 



f{t + At)^e 



(All) 



(A12) 



{,H}f^{f,H}. 



(A4) 



This second order integrator is well known as the leap-frog 
integrator. 

It is widely known that symplectic integration methods 
can achieve high-accuracy once we split Hamiltonians into 
several components. Mixed variable symplectic (MVS) 
method (Wisdom & Holman 1991; Kinoshita et al. 1991) 
splits the Hamiltonian into an unperturbed part with an 
analytic solution (i.e., Keplcrian motion when we traced 
planetary motion) and a perturbation part (i.e., mutual 
gravitational perturbation among planets) . When the sys- 
tem is nearly integrable, the Hamiltonian for the unper- 
turbed part becomes much larger than that of the per- 
turbed part, which enables the method to accomplish a 
very high accuracy in integrating the equations of motion, 
compared with conventional symplectic integrators. 

The muhiple stepsize (MSS) method (Skeel & 
Biesiadecki 1994; Duncan et al. 1998) splits a poten- 
tial into the sum of potentials of a short-range and a 
long-range forces and gives different time-steps for dif- 
ferent ranges of forces. The MSS method accomplishes 
the similar accuracy compared with the usual symplec- 
tic method applied with small time-step. The use of 
different time-steps for different interactions was pro- 
posed for the integration of molecular dynamics (Streett 
et al. 1978). GADGET-2 (Springel 2005), which employs a 
TreePM method for gravitational force calculation, adopts 
different time-steps for the long-range force derived from 
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a Particle-mesh method (Hockney & Eastwood 1981) and 
the short-range force derived from a Tree method (Barnes 
& Hut 1986). 

BRIDGE (Fujii ct al. 2007) was developed in or- 
der to solve galaxy-star cluster systems self-consistently. 
BRIDGE divides a Hamiltonian of a galaxy-star cluster 
system into a star cluster and a galaxy parts, and ap- 
plies different time-steps and integrators. In BRIDGE, 
the integration of star cluster particles is performed by 
a forth-order integration method, namely Hermit method 
(Makino & Aarseth 1992). Force calculations among star 
clusters arc performed by direct method with GRAPE 
(Sugimoto ct al. 1990). The integration of galaxy parti- 
cles is performed by the leap-frog method with the Tree 
method for the force estimation. Forces between star clus- 
ter particles and galaxy particles are also calculated by 
the Tree method with constant time-step. Therefore this 
method can deal with coevolution of collisionless and col- 
lisional systems self-consistently. 
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